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ABSTRACT 

Central gravitational image detection is very important for the study of the mass dis- 
tribution of the inner parts (~ 100 pc) of lens galaxies. However, the detection of such 
images is extremely rare and difficult. We present a 1.7-GHz High Sensitivity Array 
(HSA) observation of the double-image radio lens system B1030+074. The data are 
combined with archive VLBA and global- VLBI observations, and careful consideration 
is given to the effects of noise, CLEANing and self-calibration. An upper limit is derived 
for the strength of the central image of 180 /iJy (90% confidence level), considerably 
greater than would have been expected on the basis of a simple analysis. This gives 
a lower limit of ~ 10 3 for the ratio of the brightest image to the central image. For 
cusped models of lens mass distributions, we have made use of this non-detection to 
constrain the relation between inner power-law slope (3 of the lensing galaxy mass 
profile, and its break radius rj,. For r\, > 130 pc the power-law slope is required to 
be close to isothermal (f3 > 1.8). A flatter inner slope is allowed if a massive black 
hole is present at the centre of the lensing galaxy, but the effect of the black hole is 
small unless it is ~ 10 times more massive than that implied by the relation between 
black hole mass and stellar velocity dispersion. By comparing four epochs of VLBI 
observations, we also detected possible superluminal motion in the jet in the brighter 
A image. The B jet remains unresolved, as expected from a simple lens model of the 
system. 
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1 INTRODUCTION 

Gravitational lensing is an important method to probe the 
mass distribution in the Universe, from the largest scales 
down to sub-galactic scales. The matter in our universe 
is thought to be predominantly non-baryonic, consisting 
of a weakly-interacting form known as Cold Dark Matter 
(CDM). CDM simulations have so far been very successful in 
reproducing large-scale s tructures observed in large galaxy 
surveys su ch as the 2dF (IColless et al.ll200ll : |Peacockll2003i ) 
and SDSS (|Doroshkevich et al.ll2004l ). CDM structures form 
potential wells, and numerical simulations predict that 
these have a characteristic radial mass profile consisting of 
a broken power law with a flatter central slope, now know n 
as a NFW profile (jNavarro. Frenk. fc White] Il996l . Il997t ). 
Baryons dominate the central regions of galaxies, having 
cooled into the centres of the po t ential wells formed by 
the CDM jRees fc Ostrikerl 1 19771 : iNavarro fc Bend Il99ll : 
ICombes! |2004| ). However, the detailed physics is difficult; 
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firstly, large-scale CDM simulations reach the limit of their 
curre n t resolution on scales of ~ kpc (e.g. Ghigna et al.l 
20001; ISommer-Larsen et all 120031 : iRobertson et all 12004 



Kanget al .1120051 ) and second the baryon content requires 
the i nput of detailed micro-physics (e.g. IWhite fc Frenkl 



Il99ll : iGnedin et ai]|2004 iMaccio et aflbooel ) 

In the central regions of galaxies we have reasons to 
believe that an additional component is present, namely a 
massive black hole. At the centre of our own Galaxy there 
is str ong evidence for t he existence of a 3.7 x 10 6 Mq black 
hole (|Ghez et alj|2005h . and more massive black holes of 
~ 0.002 times the total galaxy mass have been inferred 
from kinematical observations (|Lauer et al.l Il997l ). Some 
galaxies have bla c k hol e masses as large as 1 x 10 9 Mq 
(|Kormendv et ah I 1 19961 ). Radio observations have also 
sho wn that massiv e blac k holes exist in active disk galax- 
ies (|Mivoshi et all Il995l ). The black hole mass is also 
found t o scale with the ste l lar ve l ocity dispersion of the 
bulge (iFerrarese fc Merrittl l2000l : iGebhardt etafl l200d : 
iTremaine et al. 20021 ). in a process which is probably 
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related to galaxy formation in a way that is not yet 
completely understood. 

The mass distribution in the central regions of galaxies 
can be probed by gravitational tensing, as was first pointed 
out by Wallin gton fc Naravanl (|l993l ). This is because the 
central image of a gravitational lens system passes through 
the centre of the lens galaxy potential, and its magnification 
can in principle tell us about the nature of the potential 
within ~ 100 pc from the centre of the galaxy. In general, 
the more nearly singular the potential, the more strongly 
demagnified is the central image. 

A number of theoretical and observational analyses of 
central images of galaxy lens systems have been carried out 
to date. A strict proof for the general theorem that lenses 
have odd n umbers of image s was carried out with Lorentzian 
manifolds (McKenzie 1985|), though its app licable conditions 
are arguably not universal (|Gottlieblll994h . However, most 
physical models show the existence of the central odd image 
(Section Q. Observationally, there are some cases in which 
radio emission is detected near the centre of a lens galaxy, 
for example in the radi o-loud lens systems MG J1131+0456 
dChen fc Hewitt! 1 19931 ). QSO 0957+561 dHaryanek et all 



19991) and 



1 19971 ). CLASS B20 45+265 (iFassnacht et al 
CLASS B2108+213 |McKean et al.l 120051 ). None of these 
are secure or even probable detections of a central lensed 
image, as the central component is quite likely to be due to 
radio emission from active nuclei in the lens galaxies. The 
lens system PMN J1632— 0033 has yielded the most likely 
detection of a central image, although in this case free-free 
absorption in the lensing galaxy is required to modify the 
radio spectral index of the cent ral image to make it agre e 
with that of the other images l|Winn et all 120031 . l2004bl) . 
In addition, there is one optically-selected three-image 
case, APM 08279+5255 l|lbata et al.1 Il999l ; lEgami et all 
2000; Mu noz et al.ll200ll ). which may also be caused by a 
"naked cusp" confi guration rather than a central minimum 
l|Lewis et al.ll2002T ). A likely central fifth image has recently 
been discovered in SPS S J1004+4112 in which the lens 
is a cluster of galaxies l|lnada et all 12005?). Recently the 
more systematic Extragalactic Lens VLBI Imaging Survey 
(ELVIS) project has attempted to find detections or limits 
on core images, r esulting in one n on-detection so far in 
PMN J1838-3427 jBovce et alj|2006ft . 

There are different views on the mass distribution in the 
inner regions of galaxies. Existing observations can be used 
to constrain the mass distribution of central regions of the 
lensing galaxy. Since numerical simulations predict double 
power-law models wi th central cusps rather than core radii 
l|Navarro et al.|["l997l ). and the observed light distributions 
seem also to follow cusp-like laws in the centres of nearby 
early- type galaxies (|Faber et allll997l ). these analytic forms 
are usually assumed when fitting centr al mass distributions 
|Munoz et al.ll200ll ; IRusin et al.ll2003l ) On the other hand, 
in some cases, in particular in low-surface-brightness galax- 
ies where a large da rk matter component may be present, 
Ide Blok et all (|200ll ) have argued that the central power-law 
exponent is close to flat, corresponding to a central density 
"core". Central image detection (or limits) can help us to 
understand the mass distribution. 

A thorough analysis of the available constraints from 
th e non-de t ection of t he central image h as been carried out 
bv lKeetonl (|2003l ) and iMao et "all l|200ll) have extended the 



analysis to cover the properties of central black holes. In 
principle, a non-detection of a central image implies either a 
relatively stee p inner mass slope o r a small break radius in 
the power- law. IRusin fc Mai (|200ll ). based on analysis of six 
systems, found that th e inner mass s lope index is unlikely 
to be flatter than 0.8. iKeetonl (|2003l ) discussed a range of 
models for realistic lens galaxies and found that galaxies 
which yield central image magnifications of fi ~ 0.001 are 
likely to be the most common, and that in 10-20% of cases a 
central black hole may be large enough, and positioned close 
enough to the central image, to demagnify it into invisibility. 

In this paper we cons ider the lens system B1030+074 
ijXanthopoulos et al.lll998T ). This is in many ways the ideal 
system in which to look for an odd image. First, it is very 
strong, having a primary comp onent with flux density o f 
more than 200 mjy at 5 GHz (|Xanthopoulos et al.l fl998). 
Second, being a radio source, the detection of the central 
image is made easier by the absence of strong radio emission 
from the lensing galaxy. Third, it has two observed images 
which are relatively asymmetric in flux density (~ 15 : 1). 
Such lens systems are likely to give rise to central images 
which are less strongly demagnified than fou r-image lens 
syste ms or more symmetric two- image systems (W inn et all 
2003). Our aim was to produce an image with an rms noise 
of about 10 -4 of the flux of the primary, and hence give us a 
reasonable ch ance to detect the central image based on the 
prediction of IKeetonl l|2003h . In the following Section EH 
[4] and [5] we present results from observations, a derivation 
of the upper limit to detection, results from lens modelling 
and detection of superluminal jet motions, respectively. 



2 OBSERVATIONS 

The double-image lens system B1030+074 was discovered 
during the 8.4- GHz Jodrell-VLA Astrometric Survey 
(JVA S) in 1992 (|Browne et al l Il998l ; IXanthopoulos et all 
1998). A follow-up 1.7-GHz observation with Very Long 
Baseline Array (VLBA) and Effelsberg showed the primary 
image had a flux density of 250 mjy and the secondary 
image had a flux density of 17.4 mjy with a separation 
of 1.567 arcsec at a position angle of 143°. 4. The nominal 
dynamic range of the 1 .7-GHz map was about 6000:1 
ijXanthopoulos et al.l 1200 ll ) . We have made a new L-band 
Very Long Baseline Interferometry (VLBI) observation 
with the High Sensitivity Array (HSA: VLBA, VLA, Green 
Bank Telescope and Arecibo) aimed at approaching a 
higher dynamic range in order to detect any central image 
or improve the upper limit. This 1.7-GHz dual polarisation 
observation was carried out on 2004 December 27. About 3 
hours integration were obtained with an aggregate bit rate 
of 256 Mbs -1 . In frequency this corresponds to 4 IFs each 
containing 32 channels of 0.25 MHz, using 2 bits/sample. 
The data were correlated using the VLBA correlator in 
Socorro, averaged into 2-s bins. A single field centre was 
used, centred on the mid point between image A and B. 

The data calibration was mainly carried out with AIP^|. 
Bad data were flagged, including the removal of the upper 



1 the Astronomical Image Processing System distributed by 
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two spectral channels, all baselines involving the Mauna Kea 
VLBA antenna, and a large segment of data from the Green 
Bank antenna in the first hour of observations. The VLBA 
calibration transfer, written by the correlator into the data 
file, was applied, having been corrected for an incorrect sys- 
tem temperature due to a VLA mode single dish/phased ar- 
ray recognition problem (Sjouwerman 2005, NRAO commu- 
nication). The data were re-inspected and further sporadic 
episodes of bad data were flagged. A fringe-fit was then per- 
formed using the data on the target source B 1030+074; we 
have experimented with use of a two-component model and 
a point source in the fringe fit and find that little difference 
is evident. The data after fringe fitting were checked to en- 
sure a flat response across the bandpass, and the phase and 
amplitude solutions derived during the calibration process 
were applied to the data. 

The data were averaged in time using the AIPS task 
UBAVG, which performs a baseline-dependent time averag- 
ing, giving < 1% flux reduction by time-average smearing 
over the 1600-mas extent of the source. In frequency, the 
data were then averaged into 8 channels, each with a band- 
width of 1 MHz, which gives a flux reduction of < 3% over 
the required field of view. Images were produced using the 
AIPS IMAGR program, wh ich implements the Clark v ersion 
of the CLEAN algorithm (|Hogboml [l974l ; IClarkl Il980h . The 
data were self-calibrated using the initial image as a model, 
taking CLEAN components up to the first negative. A global 
amplitude self-calibration was first applied, using a very long 
solution interval to adjust the flux scale of the antennas to 
each other. After production of a new CLEANed image, a 
phase self-calibration was then applied using a short solu- 
tion interval. After correction of phases in this way, a fur- 
ther round of self-calibration using amplitudes and phases 
was performed, with a solution interval of 15 minutes. No 
further improvement was observed in the maps by use of 
further self-calibration. Final images were produced using 
uv grid weighting close to uniform and weighting of the data 
points by the statistical weight appropriate to the size of the 
telescopes on the baseline involved. Natural grid weighting 
in principle should give better signal-to-noise by a factor of 
~ 2, but was found to give very poor results in practice. This 
is because calibration errors and CLEAN errors are more im- 
portant in high-dynamic-range mapping, we cannot expect 
the thermal noise limit to apply. In general, for producing 
the final image from the self-calibrated data we CLEANed 
the whole image, or a large part of it, using a large num- 
ber (~ 20000) of CLEAN iterations, a procedure which we 
justify in detail in Section [3] Gaussian model fitting to the 
final images was done with the AIPS task jmfit in the image 
plane. 

To get the maximum dynamic range, we also obtained 
data on B 1030+074 from previous 1.7-GHz VLBI observa- 
tions (1998 BX003, 2000 GX006, 2001 GX007), and com- 
bined them with our new observation (2004 BJ054). Data 
from the previous epochs were calibrated and imaged in the 
same way as previously described. For each epoch, CLEANed 
maps were made individually, and combined together in the 
map plane weighting by 1 /a 2 , where a is the measured rms 
noise in each image. Our combined global- VLBI + HSA map 
of B1030+074 made with a restoring beam of 7.7x4.5 mas 2 
at a position angle of —10°. 85, is shown in Fig. [T] together 



with enlarged versions of images A and B. Parameter values 
from model fitting are given in Table [2] 

Image A is resolved with an extended jet that can be 
seen in the north-east direction having a length of about 
20 mas, while image B is not resolved. A super-resolved 
image with a half-size restoring beam (i.e., 3.8x4.5 mas 2 ) 
was made, but the B jet was still not detected. The 
whole field CLEAN of the HSA map gives a noise level of 
a ~ 24.1 /LtJybeam -1 around the expected position of 
the central image (see Table [TJ. After the combination 
of the maps from four epochs, the noise level is reduced 
to a ~ 19.9 ^tjybeam -1 . The dynamic range achieved in 
the final map is more than 10000:1 (Sa/o > 10000, see 
Table [2]). 

To locate the lens galaxy and hence derive the likely 
position of the central image, we use d the HST WFPC2 
I-band image (Fig. HJ of B1030+074 (|Xanthopoulos et all 
Il998h from the HST archive to measure the position of the 
lens galaxy. We can clearly see both the lensed images and 
the lens galaxy which has an extended structure to the west. 
The separation from lens galaxy G to the B image is 127+10 
mas in position angle 149°. 7 + 4°. from north through east; 
there is no obvious detection of a radio source above the 
local noise level at the position of G. 

The flux density of both A and B components has var- 
ied (~ 10%) over the seven-year span of observations. Vari- 
ation was also seen i n a m onitoring campaign reported by 
IXanthopoulos et al.l (I2000T ) which consisted of 47 epochs of 
VLA observations spread over 240 days. Both datasets are 
not inconsistent with intrinsic variation of the source, to- 
gether with an expected time delay of ~ 100-150 days pre- 
dicted by s i mple l ens models. In the case of J1838— 3427, 
W inn et ahl (|2004aT l deduce that Galactic scintillation makes 
a major contribution to the variability; this process is less 
likely to operate in B1030+074 due to the relatively high 
(b = 52° ) Galactic latitude. 



3 DERIVATION OF THE UPPER LIMIT 

The 3-cr limit of our final HSA map is about 70 piJy beam - , 
and this level has traditionally been regarded as the effec- 
tive upper limit on the central-image flux density if there 
is no obvious detection. This is not the case, for two rea- 
sons. First, the image has undergone a CLEANing and self- 
calibration process, which may have removed flux from any 
central image. Second, the uncertainty in the position of the 
central image means that the image statistics must be con- 
sidered carefully. We consider each issue in turn. 

3.1 CLEANing and self-calibration effects 

The use of the CLEAN algorithm in the mapping procedure 
may result in elimination of f aint images of a few a, a prob- 
lem known as "clean bias" (jCondon et al.lll998r ). Simula- 
tions with artificial sources show that sources up to 6a may 
be eliminated in this way if excessively CLEANed, but that 
use of CLEAN boxes around the true source position should 
allow them to be recovered. It is also possible that the use 
of self-calibration may lead to the adjustment of the data 
in such a way as to eliminate weak sources, by partially ab- 
sorbing them into the calculated gains for each telescope. 



© 0000 RAS, MNRAS 000, 000-000 



4 M. Zhang et al. 





Figure 1. Combined global- VLBI + HSA images of B1030+074. The restoring beam size is 7.7x4.5 mas 2 at a position angle of —10°. 9. 
Labels A, B, G and S denote primary A, secondary B images, lens galaxy and the source respectively. Image A and image B are shown 
in two separate magnified figures. Contours in the maps are plotted at multiples (-1, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048) X 
3<t where a is the off-source local rms noise in the map (53 /xjybeam" 1 in A image and 26 /xjybeam - 1 in B image). The top-left figure 
shows a fitted SIE model, while the solid curves are caustics and the dashed one is the critical curve. The two boxes in the map are our 
CLEAN recovery test region and third image prediction region respectively, see details in Section [3] 



Experiment code 


Telescope 


Observing date 

(yyyy mm dd) 


Time on source 
(min) 


Measured rms noise 
(/xjy beam ) 


Thermal noise 
(/LtJy beam -1 ) 


BX003 


VLBA+Eb 


1998 06 10 


705.00 


61.8 


44.4 


GX006 


VLBA+Eb+Jb+Wb 


2000 02 12 


545.00 


58.1 


24.1 


GX007 


VLBA+VLA+EVN subset 


2001 02 09 


526.00 


39.5 


18.9 


BJ054 


HSA 


2004 12 27 


166.25 


24.1 


8.2 



Table 1. Noise levels in the central image regions of four epoch 1.7-GHz VLBI/HSA observations 





R.A. offset 
(mas) 


Dec offset 
(mas) 


Flux density 
(mjy) 


Peak intensity 
(mjy beam ) 


A 


0.0000 


0.0000 


258 ± 13 


208 ± 10 


Al 


5.613 ± 0.006 


3.158 ±0.005 


8.5 ± 0.4 


3.6 ±0.2 


A2 


20.47 ±0.03 


9.96 ±0.04 


2.6 ± 0.1 


0.73 ± 0.04 


A3 


30.0 ±0.1 


17.4 ±0.1 


1.00 ±0.05 


0.24 ±0.01 


B 


934.987 ± 0.001 


-1258.092 ±0.002 


16.9 ±0.8 


12.5 ±0.6 


G 


869.0 ± 10.0 


-1147.0 ± 10.0 







Table 2. The observational constraints from combined VLBI + HSA images. Al, A2 and A3 are jet components within component A 
(top-right figure in Fig.[TJ and G is the lensing galaxy. The position of G comes from HST I-band image (see text). 
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Figure 2. HST I-band image of B1030+074. The diffuse double- 
peak pattern at the lower part includes the B image and lens 
galaxy. Contours in the map are plotted at multiples (-1, 1, 2, 4, 
8, 16, 32, 64, 128, 256, 512) X 0.35% peak value. The two images 
and lens galaxy are labelled as A, B and G respectively. 



full flux density of the artificial sources. In Fig. [4] we plot 
the flux density inserted against the flux density recovered. 
Although there is considerable scatter in the relation, a lin- 
ear fit gives a gradient of 0.72, implying a ~ 30% reduction 
in flux density caused by the self-calibration and mapping 
process. 

One possibility is that the weighting of the individual 
antennas could be responsible for the non-recovery of weak 
sources, since the HSA is a combination of ten small an- 
tennas and three (Green Bank Telescope, the phased VLA 
and Arecibo) which have orders of magnitude greater col- 
lecting area. If the large telescopes are given their full sta- 
tistical weight, the number of constraints used to derive the 
amplitude and phase corrections applied by self-calibration 
is effectively reduced, because few baselines are weighted 
significantly and in fact the signal is dominated by three 
highly sensitive baselines. This in turn means that the so- 
lution for the amplitude and phase errors on the large tele- 
scopes is under-constrained, and allows genuine structure 
to be removed by the corrections. To test this possibility, 
we weighted down the large telescopes, up to the extreme 
case of treating all telescopes with equal weight, at a cost of 
an increase in the thermal noise. However, we did not find 
much improvement of the recovery of the added sources, and 
in fact most added sources suffered worse flux loss. 



In order to assess this problem, we have used the 
AlPS task uvmod to add artificial point sources with a 
range of fluxes to the calibrated uv data, and determined 
whether they are recovered during passage through our 
self-calibration and mapping procedure. 18 artificial point 
sources were used simultaneously, arranged in a ring of 
diameter 200 mas (see Fig. [3} and with flux densities 
ranging from 60 /iJy to 1 mjy. This also allowed us the 
opportunity to adjust the procedure in order to recover 
the weakest possible source. The philosophy throughout 
the process has been to apply the same procedure to the 
artificial images as to the possible central image. 

For reasons which we have been unable to determine, 
the use of relatively tight CLEAN boxes around components 
A and B gives poor noise characteristics and a high rms in 
the rest of the image, including the position of the possible 
central image. On the other hand, CLEANing the whole image 
gives the low noise (24 ^Jybeam -1 ) seen in Fig. [T]but rela- 
tively poor recovery of the artificial images. An alternative 
boxing scheme is to CLEAN the majority of the image, leaving 
a gap around the ring of artificial sources and the position of 
the central image, but in addition including a small box of 
10 x 10 mas 2 around the expected position of each individual 
artificial image and of the central image. This gives good re- 
sults in that artificial sources of 100 ^iJy can be recovered; 
because the artificial sources and any potential central im- 
age are treated identically, we can be confident that any real 
central image at this flux density level could be detected. 

The results are shown in Figs. [3] and [4] The best recov- 
ery is from our HSA data, with recovery of 100 /xJy sources. 
However, in the previous three VLBI observations, stronger 
sources than this (up to 200 fJ,3y) were in many cases not re- 
covered. In both cases, the recovery of artificial test sources 
is only achieved at between 5 to 10 times the expected ther- 
mal noise. In addition, our procedure does not recover the 



3.2 Statistics of non-detections 

We will take a simple and conservative approach to assessing 
the detectability of a third image in our data. We would like 
to be able to give a value for probability that we would 
have detected a third image of flux density ^ 5* in our data, 
if it were there. Unfortunately we cannot do this. We can 
assess the probability of detecting a source of flux density 
S by doing simulations but not the cumulative probability 
of detecting an image ^ S. This is because we do not know 
the expected distribution of the third image flux densities 
over which we should integrate. However, one thing we can 
be confident of is that the detectability of images increases 
with increasing flux density. This means that, if we assess 
the probability P* of detecting images of flux density S by 
adding artificial sources of this flux density to our observed 
data, and see how many we recover, the result can give us 
a lower limit on the desired cumulative probability. We will 
adopt P* =0.9 as our conservative limit on the probability 
of detecting an image, if it had a strength of > 5*. We now 
describe how we search the dataset before we return to the 
probability calculation. 

Because of uncertainties in the galaxy position, and be- 
cause the third-image position depends on the adoption of a 
model, we have to consider the pixel values in the final map 
over a search box rather than at a single position. 

Defining the search box depends on the adoption of a 
model. We defer detailed discussion of this model to the 
next section, but note that its two principal free parameters 
are the break radius rt and inner power-law slope (5 of the 
galaxy mass distribution. Fig. [5] shows the parameter space 
explored to find possible positions of the central image, and 
the predicted third image positions are shown in Fig-EI From 
Fig. [6] we found that all the possible third image positions 
are enclosed in a club-shaped pattern along the major axis 
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Figure 3. Artificial source 'ring' recovery test. The solid rectangles are used as a guide to the source positions. The positions of our 18 
added models in the uv data form a circle with 20° separation from each other. The upper half-circle includes 10 sources from 1000 ^Jy 
to 100 ^Jy in steps of 100 £tJy, and the lower half-circle includes 8 sources from 200 /ijy to 60 /^Jy in steps of 20 /^Jy. For comparison, 
all maps are set to linear scale with a flux density range of from 0.03 mJy to 1 mJy. 
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Figure 4. Flux loss in the HSA map from our mapping pro- 
cedure, expressed in a plot of the flux density of an artificially 
added source against the recovered flux density. The dashed line 
is a linear fit to the data which has a slope of 0.72. The dotted 
diagonal line has a slope of 1. 



of the elliptical mass distribution. The area of the club shape 
is about 840 mas 2 , i.e., of about 24 beams. 

Using 10x10 mas 2 CLEAN boxes in this area for imag- 
ing, as described in the previous subsection, gives a maxi- 
mum flux density within this region of 42 /ijy. This value is 
consistent with being a noise spike, and in fact if the noise 



were Gaussian (which it is not) one would expect several 
such values within the 840-mas 2 search box. 

We now return to the probability calculation. We re- 
cover no source stronger than 42 fiJy from the search box. 
What is the flux density, S, that an image would have to 
have for P* = 0.9? To find this we have experimented with 
artificial sources of different flux densities and find that us- 
ing sources of 180 ^Jy, 90% of them are recovered with flux 
densities greater than 42 /ijy. Therefore our conclusion is 
that the probability of detecting an image using our search 
technique is 90%, provided its flux density is ^ 180 /ijy. 
A more detailed interpretation of this upper limit is given 
in Appendix 17X1 



4 LENS MODELLING 

The mass distribution of galaxies has been a subject of de- 
bate in the last two decades. On the scale of the outer images 
in lens systems (normally a few kiloparsecs) an isother mal 
profile is usually a good fit (e.g. iKoopmans et al]|2006l ). In 
the past, power- law models with a near ly isotherma l profile 
and a central flat core were used fe.g. | Lauerlll985T ). How- 
ever, numeric al simulations ([Navarro et al.lll997l ) and HST 
observations dFaber et al.ll 19971 ) have been used to argue for 
a more general description in which galaxies have a two- 
power-law form with a central cusp. Normally, for a weak 
cusped mode l, the odd number theorem applies and a third 
image exists l|Evans fc Wilkins"o"nlll998l '). 

We adopt a typical spherical density profile of a cusped 
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Figure 5. Grid-search sampling in /3 — r b plane. The '+' symbol 
denotes the grid points we sampled to determine the predicted 
central image position. Above the dotted line, the central image 
is predicted at a level which would be obvious from inspection of 
the radio map. 



model fe.g. lMunoz et al.ll200ll ): 

p(r) 



Po 



{r/r b )0[l + (r/ri,) 2 ]("-«/ 2 



(1) 



where p is the mass density, po is the scale factor, r is the 
radius, r\ is the outer slope, is the inner slope, and r b is the 
break radius. This two power-law distribution degenerates 
to single power-law distributions r~ v or r _f3 respectively 
when r 3> r b or r -C r b . In the limit — » 0, we have a 
cored model, whereas in the limit = rj = 2.0, we recover 
a singular isothermal model (SIE). More generally, for the 
ellipsoidal models, the r should be repl aced by triaxial c o- 
ordinates (r 2 = x 2 /a 2 + y 2 /b 2 + z 2 /c 2 ) (|Chae et alJll998h . 

Numerous authors have pointed out that regions of the 
r b -0 parameter space can be ruled out by the non-detection 
of central images in gravitational lens systems. In all cases a 
power-law model with an isothermal slope or steeper causes 
inf inite demagnification of the c entral image. Early studies 
bv lWallington fc Naravanl (| 19931 ) examined softened power- 
law models with a core; in t his case an upper limit on the 
core radius could be derived. lEvans fc Hunter! l|2002h devel- 
oped a general formalism for cored profiles and found that 
for observed systems with missing central images, the typi- 
cal core radius of the lensing galaxy must be <300 pc; they 
concluded that the mass profile must be nearl y cusped and 
that t he cusp must be isothermal or stronger. iRusin fc Mai 
|200ll ) used the non-detection of central imag es in a num- 
ber of Cosmic Lens All-S ky Survey (CLASS) dMvers et al.l 
120031 ; [Browne et al.ll2003l ) gravitational lens systems to infer 
limits on single power-law models (in which = 77 and r b 
is therefore irrelevant) and for B1030+074 found that the 
then existing limit on the central image imp li ed a central 
power law slo pe > 1.91.lMunoz et all (|200ll ) , IWinn et all 
(|2003l . l2004ah and lBovce et al l (|200fj ) explicitly studied the 
allowed region of the r b -0 plane for various lens systems, 
assuming an isothermal outer slope (?7=2). In each case the 
non-detection of a central image allows the exclusion of a 
region of this plane with a relatively flat inner slope (low 
0) and large break radius (high r b ), and this restriction be- 



comes more powerful as the limit on the ratio between the 
central image and the brightest i mage becomes mor e strin- 
gent. In the case of a detection (|Winn et al. 2004b) an al- 
lowed locus on this plane can be inferred. 

The number of degrees of freedom in a simple cusp 
model is negative, since the free parameters (dimensionless 
surface density kq, ellipticity e, position angle 8, r b , 0, 
source and galax y (x, y) positions) o utnumber t he measured 
param eters (e.g. iBovce et all 120061 ) . Following IWinn et al.l 
|2003t ) and IBovce et alj (|2006t ). we therefore investigated 
the behaviour of cusped models for B1030+074 with differ- 
ent values of inner power-law slope, 0, and break radius, 
r b , together with a fixed isothermal outer power-law slope 
77 = 2. This was repeated for different flux ratios between 
the primary image A and the third image C ( see Fi g. 0. 
We us ed th e fast Fourie r expansion method of IChae et al.l 
(|l998h and IChael (|2002T ) for this purpose. The higher the 
observational limit on the magnification ratio, the larger 
the non-detection region above the solid lines. For our value 
of Sa/Sc ~ 1000 (from our derived upper limit on So) and 
r b > 130 pc, we obtain an inner power-law slope greater 
than 1.8. 

One possibility is that a massive central bl ack hole in 
the lens galaxy could suppress the centra l image (|Mao et al.l 
l200ll ; iKeetonfeoOll . 12003 ) ■ iKeetonl (|2003l ) analysed the likely 
suppression of central images for different black hole masses 
and concluded that the average degree of suppression is quite 
small; the exception is in those cases where the cen tral image 
is extremely faint in any case. IRusin fc Mai l|200ll ) calculate 
explicitly the case of B1030+074, using a single power-law 
model (0 = 77) together with a plausible black hole mass, 
and conclude that the limit on the power-law slope needs to 
be relaxed from > 1.91 to > 1.83. 

A plausible black hole mass can be estimated by 
using the empirical correlation between the central 
black hole mass and the galact ic bulge velocity d isper- 
sion dFerrarese fc Merrittl 120001 : ICebhardt et all l200fj| ; 
iTremaine et al.ll2002l )~ 

log(M BH /M e ) = a + 0\og(a/a o ), (2) 

where = 4.02 ± 0.32, a = 8.13 ± 0.06, a = 200 kms" 1 
are the fitted values from ITremaine et~ al. (2002). Since the 
galactic mass distribution is very close to isothermal, we can 
use the relation 



A6/2 = in(a 2 /c 2 )(D d3 /D 3 ) 



(•3) 



to estimate the velocity dispersion of the lens galaxy, where 
A8 is the angular separation between image A and B, Dd s 
and D„ are the angular diameter distances between lens and 
source, and source and observer, respectively. For the case of 
B1030+074, A9 ~ 1.567 arcsec, which gives a velocity dis- 
persion of 233 kms -1 , and hence a central black hole mass 
Mbh of about 2.5 x 10 s Mgjf]. Applying the cusped model 
together with this black hole mass, we can get the modified 
relation between break radius r b and inner power-law slope 
0, which is shown in Fig. [7] In the limit of large break radius 
(r b — > 00), the cusped model app roaches a single pow er-law 
model similar to the simulation o f R usin fc MaT(|200ll ) previ- 
ously described. In our model, the asymptotic value of inner 

2 Here and elsewhere in this paper we use a flat ACDM cosmolog- 
ical model with S7 m = 0.3, Ha = 0.7 and Hq = 72 kms _1 Mpc _1 . 
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Figure 6. Predicted third image positions from cusped models. The solid circle is the galaxy position with lcr error. The dot-dashed line 
denotes the major axis of the elliptical mass distribution. In the left figure, the '+' symbols denote the predicted third image positions 
in the image plane from each fit using different values of rj, and (3. The club-shaped dashed line pattern encloses all the predicted third 
image positions in our test and has an area of 840 mas 2 ; the top and right panels show the \ 2 distributions of each fit, as a function of 
ARA and ADEC respectively. In the right figure, we show the HSA map of the predicted third image region overlapped with the club 
pattern, and the rms uncertainty of the galaxy position as a circle. In these models, the galaxy position has been allowed to vary, but 
with a Gaussian penalty function based on our knowledge of its position from the HST image and the rms uncertainty. 



power-law slope f3 ~ 1.95, and /3 ~ 1.91 with a 2.5 x 10 s M 
central black hole in the lensing galaxy. If we cut off the 
break radius at the Einstein radius of an isothermal model, 
we get an inner power-law slope j3 ~ 1.93, and /3 ~ 1.88 with 
a 2.5 x 10 s Mq central black hole. Obviously, a central mas- 
sive black hole affects the rj, — j3 curve more when the flux 
ratio Sa/Sc is bigger. The 2.5 x 10 8 M Q black hole derived 
from the Mbh — a relation does not soften the cusped mass 
profile significantly. We would need to add a black hole of 
10 times greater mass in order for the cusped mass profile 
to be significantly shallower for the same Sa/Sc constraint. 
This result agrees well with Keeton's conclusion from sim- 
ulations of star+halo CDM mass models, that black holes 
would have to lie off the Mbh ~ a relation by at least a 
factor of 10 in mass to significantly affect the central images 
|Keetonll200ll ). 



5 SUPERLUMINAL JET? 

The jet within the A component is potentially useful 
astrophysically. First, because of the stretching induced 
by lensing magnification, any superluminal motion is 
easier to observe. Superlumin al motion was predicted 
four decades ago (|Reesl 1 19661 ). and many superluminal 
jets have been discovered since then l|Vermeulen fc Cohenl 
1 19941 ; iBritzen et al.l Il999i ) . These superluminal motions are 
typically a few times the speed of light, but would appear 
larger in the image plane because of magnification; in 
extreme high-magnification cases speeds up to 100 c would 
in principle be observable l|Hogg fc Blandfordl [l994) . The 
nearly seven-year span of our observations of B 1030+074 
gives us a good opportunity to investigate the existence of 
any superluminal movement of the jet in a lens system (see 
Fig. [8}. Second, any components found in both A and B 



give extra information to constrain the mass model. The A 
image in the B1030+074 lens system has a very significant 
jet while the jet is not resolved in the B image. 

Maps of image A at all four epochs, made using the same 
CLEAN restoring beam, are shown in Fig. [8] At each epoch 
the jet exhibits a number of knots. Although there are always 
ambiguities in identifying such jet components from epoch 
to epoch, the most conservative correspondence is shown by 
the dotted lines in Fig. [8] There is a suggestion of outward 
motion of the outermost jet knot 3; the fitted slope implies a 
proper motion of about 0.34T0.20 mas yr _1 , corresponding 
to a velocity of the jet component in the image plane of 
about 7.2 ± 4.2 c. The fitted magnification of the A image 
for an SIE model is about /j, ~ 3.0 in the direction along the 
jet, so in the source plane the jet proper motion is about 
3. Oil. 7 c. From top to bottom, the other two jet components 
give proper motions in the source plane of 3.7 ± 2.8 c and 
1.1 ± 0.9 c, respectively. The determination of the true jet 
proper motion relies on the magnification factor which is 
model dependent. The cusp model can give a range of A 
image magnification from 3.0 to 5.0. With the Sa/Sc ~ 
1000 non-detection constraint the fitted cusped model gives 
a typical A image magnification of about 2.8, and 3.1 with a 
2.5 x 10 8 Mq black hole estimated from the Mbh—o relation. 
A larger magnification of image A could be achieved if we 
add a bigger black hole mass to suppress the central image. 

The jet in the B image is not resolved. Using the com- 
ponents of the observed A jet from the VLBI/HSA image 
together with the fitted lens model, we simulated the struc- 
ture of the B jet using a point spread function (PSF) similar 
to the restoring beam in our VLBI images (Fig. [5| . The sim- 
ulation confirms that we would not expect to resolve the B 
jet structure using our current observations. 
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Figure 7. Third image detection regions in break radius - in- 
ner power-law slope diagram. From the top to the middle and 
the bottom, the three figures are calculated with primary /central 
image magnification ratios of 10 3 , 10 4 and 10 5 respectively. The 
solid lines are the results calculated by the cusped model, the 
dashed lines are the results calculated by the cusped model with 
a 2.5 X 10 8 Mq black hole, the dot-dashed lines are the results 
calculated by the cusped model with a 2.5 X 10 9 Mq black hole. 
Re is the major axis from a fitted SIE model. 
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Figure 8. Jet proper motion of B1030+074 image A. The top 
figure shows the most likely correspondence of outgoing indi- 
vidual jet components. All images are restored with the same 
restoring beam and aligned assuming a stationary core position 
where the vertical scale is proportional to the time intervals be- 
tween the different observations. The bottom figure shows the 
distances between jet components and the core at a position an- 
gle of 60°. 1. The labelled numbers show the correspondence of 
jet components. The three dashed lines indicate the results from 
weighted fits to proper motions of the jet components (2a, 2b, 
3), corresponding to the three dashed outlines in the top fig- 
ure. Fitted slopes are 0.34 ± 0.20 mas yr — 1 for component 3 



be 



0.07), 0.42 ± 0.32 mas yr" 1 for component 2a (x 2 = 0.7) 



and 0.12 ± 0.11 mas yr _1 (x 2 = 7.6) for component 2b. 
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Figure 9. Simulated images from the fitted SIE model. They are 
convolved with a restoring beam of 7.7 X 4.5 mas 2 at a position 
angle of —10°. 85. The contour levels here are the same as Fig. [l] 
The jet structure in image A is detectable while image B still 
looks point-like. 



6 CONCLUSIONS 

We have presented a new HSA 1.7-GHz observation of the 
double-image lens system B1030+074 which was aimed at 
finding a third image to help to understand the central 
mass distribution of the lens galaxy. Careful consideration 
of observational and statistical effects shows that the 90% 
confidence level of the non-detection in such cases is about 7 
times the rms noise (in this case corresponding to 180 /iJy). 
This gives a minimum flux ratio of about 10 3 between the 
brightest image and the central i mage. By comparing with 
predictions from a cusped model (iMufioz et alj |20011. this 
non-detection constrains the break radius - inner power-law 
(rt — /3) relation, which excludes most of the rt — region 
depicted in Fig. [S] and implies a steep inner power-law 
slope close to isotherm al. Following the Mbb — & relation 
l|Tremaine et al.l |2002| ). we predicted a central black hole 
mass of about 2.5 x 10 s Mq. However, including this black 
hole mass into the cusped model did not affect the inner 
power-law slope very much, and would only do so if the 
black hole mass were increased by a factor 10. The result is 
consistent with sim ulations from CDM mass distributions 
of the lens galaxies (|Keetonl l2003), 

Four-epoch observations gave a possible detection of su- 
perluminal motion in the jet. The determination of the true 
superluminal velocity also depends on the mass model, but 
with the fitted cusped model, the superluminal motion of 
the jet is about 3 c in the source plane, although subject to 
large errors at this stage. 

The current 1.7-GHz HSA observation cannot resolve 
the jet in the B image even with a super-resolved restoring 
beam. We simulated the B jet from the fitted mass model 
and convolved it with the restoring beam. This simulated B 
image still looks like a point source, which is consistent with 
the observation. A 5-GHz observation may possibly resolve 



it if the spectral index is not so steep as to make it below 
the surface brightness detection limit. 

In the next few years, interferometer arrays will 
come online which should make the detection of central 
images in gravitational lens systems routine. In particular, 
the e-MERLIN array will have an L-band sensitivity of 
4.0 /^Jybeam -1 with 12 hours on-source integration. With a 
resolution of 50 mas at 5 GHz, this instrument is especially 
well-matched to the problem, since the simpler image struc- 
tures are more easily mapped with high dynamic range, 
and the region to be searched to find the central image is at 
most a few res olution eleme nts. If current theoretical work 
is correct (e.g. Kccton 20031) this should allow detection of 
central images, and therefore greatly improved constraints 
on mass models, in many known radio lens systems. 
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APPENDIX A: ABOUT THE UPPER LIMIT 

There are two schools of inferential statistics: frequentist and 
Bayesian. Each gives different interpretation of probability 
as well as confidence level. We will investigate two different 
views on the upper limit. 

In the following discussion, we use the following nota- 
tions: the conditional probability P(S\D) is for "getting the 
source, given the data", and P(D\S) for "getting the data, 
given the source"; P{D) and P(S) are for data and source 
probability distributions respectively. 

Our specification of the upper limit with desired 
confidence le vel is directly from the standard Neyman 
construction l|Nevmanlll937h . To set an upper limit on our 
non-detection, we simulated the input with artificial sources 
and analysed the output from our mapping procedure. 

Our test is based on several factors: 

o The noise in a CLEANed map is not an objective statis- 
tic. Though we can manage to CLEAN it down to thermal 
noise, it is pointless if the expected weak source is CLEANed 
out. 

o The recovery of an injected artificial source is station- 
ary and related to its input position and strength, i.e., if we 
repeat trials with the same input position and strength, the 
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recovery of the artificial source through the same mapping 
procedure is same, either false dismissal or not. 

o The higher the input strength, the better recovery we 
get, while the lower the input strength the worse recovery. 
For a source above 250 fiJy, the recovery is almost certain 
though it will definitely suffer a lossy CLEANing; below f 00 
/j,3y, the recovery is hardly distinguishable from noise. 

o The background noise of a CLEANed map is non- 
Gaussian. In sidelobe contaminated region, it is rather 
obvious. To be unbiased, the recovery test should be in the 
predicted central image region. We must compromise since 
we need equal treatment (clean them simultaneously) for 
both the central image and the artificial sources. We chose 
the region close to the central image region to minimise 
this position dependency effect. To avoid confusion with 
false alarms, we should avoid inject artificial sources on 
noise spikes. This is easy to check without artificial source 
injection, and for the small sample test the false alarms are 
negligible. For example, for Gaussian noise, if you inspect 
20000 beams, you will expect one 4a noise spike. 

o There is no strict constraint on the expected central 
image flux density distribution unless it violates the flux 
ratios allowed by the cusped lens model, i.e., the central 
image should be less bright than the saddle point image B 
(Sc < Sb)- Of course, we need to be objective since we can 
always resort to more complicated models if something odd 
happens. 

o A null hypothesis. Since we found no prominent image 
in the central image region and there is no way to tell it apart 
from noise spikes even if it is present, we claim non-detection 
of the central image in our observation. We wish to derive 
an upper limit on the source strength. 

o A suboptimal strategy. The optimal statistic could 
be derived according to the Neyman- Pears on criterion if we 
know the accurate mathematical model of the noise of CLEAN 
l| Allen et al.ll2002T ). However, in the absence of any nois e 
analysis of the CLEAN procedure itself (jCornwell et al.l l 1999). 
we conservatively choose the maximum statistic (noise peak) 
from the central image region in the recovery test as our sub- 
optimal statistic, accompanied by our null hypothesis. This 
statistic is taken as the counting threshold in our recovery 
test. 

Since we found from our "ring test" that false dis- 
missals happened with artificial sources less than 200 fiJy, 
we injected 20 artificial sources of 180 /J,3y with varying 
positions close to the central image region and did the 
recovery test. Most of them are recovered with flux losses 
of ~ 20-30%, but two of them are recovered with maximum 
fluxes lower than the 42 /xJy noise peak we obtained from 
the central image region. 

(i) Frequentist upper limit: 

This is an "objective" approach which is mainly based 
on the counting frequency in our test. In fact, the standard 
Neyman construction is a frequentist method. This method 
has been extensively used in weak detection expe riments 
in particle physics and gravitational wa ve physics (|Barlowl 
l2003l ; ICranmerll2003l ; lAbbott et al.ll2004h . 

Let D be the detection statistic, and Do the suboptimal 
statistic of a targeted search in our test. If we inject artificial 
sources greater than or equal to So in the recovery test and 
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Figure Al. Central image CLEANing and artificial source recov- 
ery counting. The central image region is covered by 10 X 10 mas 2 
CLEAN boxes. The statistics Do obtained from imaging is used as 
the counting threshold. The right panel shows the possible false 
dismissals with artificial source injection and false alarms without 
artificial source injection in the real test. 



we get a fraction C of trials having residuals above Do, then 
we have a frequentist upper limit So(C) on the strength of 
the targeted source, with confidence level C. This frequentist 
confidence level should be calculated from the integrals of 
the probability distribution function (PDF) p(D|So): 



C(So)=P(D^Do|S ) 



p(D|S* )dD 



p(D\S )dD 



(Al) 



In our test, we can count how many values of D are 
greater than or equal to Do and divide it by the total number 
of D values. This is illustrated in Fig. IA1I Our null and 
alternative hypotheses are 

H : D(S) < Do, 

Hi : D(S) > Do, (A2) 
and the false alarm and false dismissal probabilities are 



a(S) = / p(a;|Ho)da;, 
JxeH! 



13(3) 



p(x\Hi)dx. 



(A3) 



We denote the detection probability j(S) = 1 — a(S) — f3(S), 
then the confidence level will be the detection probability at 
a source strength So in this case 

C(So) = 7(5*0) = 1 - o(5 ) - P(So) m 1 - P(S ). (A4) 

Since we have 2 of 20 (remember the predicted central image 
region is of about 20 beams) 180 /j,Jy artificial sources cov- 
ered above 42 pjy in our test, we will therefore say that if 
the source to be detected is at 180 /iJy level, it should be de- 
tected in our map with a flux over 42 pjy, with a confidence 
level of (20 - 2)/20 = 90%. 

We have found from the recovery test that weak sources 
are more likely to be CLEANed out, while strong sources more 
unlikely, i.e., 

P(D > D |S < So) < P(D >Do\S = So), 

P(D > D |S > So) > P{D > D \S = S ), (A5) 
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then according to the sum rule 

P{D > D \S < So) < P(D >D \S = So), 
P{D > D |S > So) > P(£> > D |S = So). 



(A6) 



So if there is a true source of 180 /xJy of higher, the the 
probability of getting an observation as small as we obtain 
or smaller is only 10% or less, i.e., to argue that there is 
a strong source which happens to have a low fluctuation 
is very implausible. If we deny the existence of a ^180 /iJy 
source in this data, we will be right at least 90% of the time. 

One should be aware that at different targeted source 
strength, we could get different confidence levels. Moreover, 
our recovery test could be biased because of the sampling 
size and position dependency. 

(ii) Bayesian (conditionalist) upper limit: 

The Bayesian approach has been consistently used in all 
statistical studies, especially for hypothesis testing of the- 
oretical predictions provided we know the prior probabil- 
ity distributions of certain physical quantities. This "subjec- 
tive" appro ach (for informative object ive Bayesian interpre- 
tations, see ljavnes fc BretthorstlboOSt ) which always needs a 
prior knowledge of a certain probability distribution, might 
not always be practical. The Bayesian confidence level ( "de- 
gree of belief") of the upper limit is calculated from the 
posterior probability P(S ^ Sq\D ^ -Do)- We can invert the 
probability P(D\S) according to Bayes' theorem 



P(Ss:S \D > Do) = 



P(D ^ D \S ^ S )P(S ^ So) 



So 



P(D > D ) 
p(S\D > D )dS 

poo 

/ p(S\D > D )dS 
Jo 

f ° P {D > D \S)p{S)dS 
Jo 

/ p(D > D \S)p(S)dS 
Jo 



(A7) 



Now to calculate the posterior probability P(S ^ So\D ^ 
-Do), we need knowledge about the probability distributions 
P(S < So) and P(D ^ D ). However, the probability 
distribution of real source strength p(S) is not known to 
us. Naturally we could assume a uniform (improper) prior 
here p(S) = constant. Even though we do not have much 
information about the data bias due to CLEANing, we can 
integrate the total probability of P(D) decomposed by the 
"cause" probability p(S), see Eg. IA7I After normalisation, 
the probability P(S So|-D ^ -Do) calculation becomes the 
probability P(D ^ Do|S ^ So) calculation, i.e., counting 
the number of recoveries in the recovery test. Indeed, 
this uniformly most powerful (UMP) test does exist if we 
assume a bound S range, i.e., the integration p(-D|S) over 
(0, oo ) converges (for m ore descript ion of tests in absence of 
a prior distribution, see ISelinlll966l ). However, the Bayesian 
confidence level from the UMP test will be very different 
from the frequentist one. Furthermore, the truth is that 
the unknown prior source strength distribution we assumed 
will affect the confidence level determination significantly. 
Since we have explored from our recovery test that at 
higher flux density level the recovery probability is higher, 




ISO 200 300 
Flux density (^Jy) 

Figure A2. Bayesian upper limit. The solid curves present 
the prior p(S), likelihood p(D\S) and the posterior probability 
p(S\D); the dashed curves present another case of different prior 
and consequently different posterior probability of which the con- 
fidence level CL increases. 



at lower flux density level the recovery probability is lower, 
we assume a monotonic continuous PDF of p(D\S) which is 
illustrated in Fig. |A2] We can see from Fig. |A2] and Eq. [A7l 
that the probability density p(S\D) is totally dependent 
on the profile of prior source PDF p(S). The unnormalised 
inverse probability density p(S\D) is the product of p(D\S) 
and p(S). So the Bayesian confidence level of the upper 
limit will be the fraction of the shaded part in the whole 
region below the curve, i.e., the cumulative probability 
P(S ^ So|-D). In the illustrated case, the confidence level is 
definitely less than 90%. Obviously, we can manipulate the 
shape of p(S) to make the confidence level arbitrary. 

In brief, if we want to set up an "objective" upper limit 
from an experiment, without assuming a Bayesian prior, we 
should subscribe to the frequentist point of view. The ap- 
proach taken in our paper is indeed a frequentist approach. 
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